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1.0  Introduction 


Considerable  research  effort  has  been  expended  in  the  general  area  of  time-frequency 
signal  analysis  tools,  and  a  wealth  of  published  material  exists  as  a  result  of  these  efforts.  To 
date,  this  research  can  be  characterized  as  a  multiplicity  of  interrelated  approaches  to  address  the 
need  for  such  tools.  Most  if  not  all  of  the  approaches  use  some  combination  of  three  basic  signal 
processing  techniques:  i)  filtering,  ii)  demodulation,  and  iii)  decomposition.  The  research 
reported  herein  falls  into  the  latter  category  of  signal  decomposition. 

Included  among  the  approaches  that  became  widely  used  by  practitioners  is  the  band-pass 
filtering  and  envelope  demodulation  approach.  This  approach  became  particularly  popular  in 
hardware  instruments  provided  by  test  equipment  manufacturers,  due  to  its  simplicity  and 
effectiveness.  The  more  popular  example  is  that  of  the  spectrum  analyzer  used  by 
communications  engineers.  These  devices  essentially  consist  of  a  fixed  frequency  narrow-band 
band-pass  filter  that  precedes  an  envelope/amplitude  demodulator.  The  input  signal  is  translated 
in  frequency  using  nonlinear  mixing  with  a  swept  local  oscillator  prior  to  band-pass  filtering. 
Rapidly  repeated  sweeps  thereby  allow  for  an  “instantaneous”  determination  of  the  frequency 
content  across  the  frequency  band  of  interest. 

From  a  research  perspective,  a  sub-category  of  approaches  has  become  very  popular  and 
can  be  characterized  by  their  use  of  instantaneous  frequency  concepts.  (For  a  single  carrier 
component  modulated  in  either  phase  or  frequency,  measurement  of  the  instantaneous  frequency 
constitutes  FM  demodulation  [Noga].)  These  approaches  are  generalized  in  the  sense  that  for 
signals  modeled  as  consisting  of  multiple  modulated  carrier  components  existing  simultaneously, 
the  goal  is  to  demodulate  each  of  these  components.  The  instantaneous  frequency  of  each 
component  is  of  particular  interest  as  a  time-frequency  (or  more  correctly,  time-instantaneous 
frequency)  representation  of  the  input  signal.  This  is  an  important  distinction  to  make  regarding 
the  various  signal  analysis  tools  that  exist.  Researchers  that  employ  techniques  from  this 
instantaneous-frequency-based  sub-category  are  often  interested  in  using  the  resulting  time- 
frequency  representation  to  detennine  how  the  signal  was  created.  For  example,  analyzing  a 
music  signal  consisting  of  the  recording  of  a  flute  may  allow  for  automatically  determining  the 
sequence  of  notes  that  were  played  to  create  the  music.  In  this  form  of  signal  analysis,  the  true 
spectral  content  is  not  necessarily  of  interest.  As  another  example,  consider  similar  analysis  of  a 
frequency-shift-keyed  (FSK)  communication  signal.  For  communication  purposes,  determining 
the  shift  frequencies  as  a  function  of  time  is  analogous  to  detennining  the  notes  as  a  function  of 
time  in  the  previous  music  signal  example.  However,  communication  designers  have  to  work 
under  various  constraints,  including  constraints  on  bandwidth  and  spuriously  generated  signals. 
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Therefore,  although  demodulation  is  of  interest  at  the  intended  receiver  of  the  FSK 
signal,  spectral  content  outside  the  location  of  the  instantaneous  frequencies  is  also  of  interest  to 
prevent  interference  to  other  users  in  adjacent  channels.  For  the  music  signal,  an  example  of 
where  the  true  spectral  content  is  of  interest  would  be  in  analyzing  tonal  quality  of  the  flute.  It 
can  be  expected  that  variations  in  tonal  quality  would  be  reflected  not  only  in  the  instantaneous 
frequencies,  but  also  in  the  spectral  regions  outside  these  frequencies. 

Popular  approaches  to  signal  decomposition  can  also  lead  to  time-frequency 
representations.  Example  methods  of  signal  decomposition  include  wavelet  decomposition,  and 
sinusoidal  basis  function  decomposition  as  when  using  the  discrete  Fourier  transform  (DFT)  or 
the  fast  Fourier  transform  (FFT).  Of  particular  interest  in  this  report  is  the  FFT-based  method 
referred  to  as  the  spectrogram.  Due  to  its  computational  simplicity  and  relative  output  quality, 
the  spectrogram  has  become  a  highly  favored  tool  for  spectral  analyses.  The  spectrogram  of  a 
signal  is  calculated  starting  with  the  selection  of  a  signal  segment  length,  L.  The  sequential 
values  resulting  from  uniform  sampling  of  the  signal  (i.e.,  the  sampled  signal)  are  separated  into 
consecutively  indexed  segments.  Each  segment  contains  a  subset  of  sequential  samples  from  the 
original  sequence.  The  subsets  are  chosen  to  start  at  sample  indices  that  increase  as  the  segment 
index  increases,  and  are  often  overlapping.  (In  this  report,  a  50%  overlap  is  used  and  segment 
indices  start  at  one  and  increase  by  0.5  for  axis  labeling  purposes.)  With  the  selection  of 
segment  length  and  the  separation  into  segments,  an  FFT  is  calculated  for  each  segment.  The 
resulting  magnitude  values  for  each  frequency  bin  from  consecutive  segments  are  plotted  versus 
segment  index  to  form  the  spectrogram. 

The  published  results  in  the  above  areas  are  both  numerous  and  high  in  quality.  General 
reference  is  made  to  [Cohen],  [Oppenheim]  and  [Suter]  for  background  in  some  of  the  various 
existing  approaches.  More  directly  related  background  for  the  research  in  this  report  can  be 
found  in  [Kay]. 

1.1  A  Short  Segment  Fourier  Transform 

For  direct  comparison  to  the  spectrogram,  the  same  segmentation  process  previously 
described  is  used  to  generate  the  new  time-frequency  representation  developed  in  this  research. 
Because  the  new  technique  also  represents  an  estimation  of  the  true  Fourier  transfonn  of  each 
segment,  it  is  by  definition  continuous  in  the  frequency  variable.  To  further  allow  for  direct 
comparisons  and  for  practical  purposes,  the  estimated  Fourier  transform  of  each  segment  is 
sampled  in  frequency  at  uniform  spacing  over  the  [0, 2n)  frequency  interval. 
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The  intent  of  the  new  approach  is  to  overcome  a  characteristic  problem  associated  with 
the  FFT-based  spectrogram.  Specifically,  as  segment  lengths  decrease  for  more  time  resolution, 
frequency  resolution  degrades  in  the  spectrogram.  An  approach  similar  to  that  taken  here  is 
presented  in  [Kay].  Important  differences  will  be  summarized  in  the  concluding  section  of  this 
report. 
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2.0  Reduced  Bandwidth  Representation 


To  develop  a  reduced  bandwidth  (Fourier)  representation  of  any  given  signal  segment, 
we  first  consider  the  response,  z/(z  —  a),  for  an  example  pole,  a  =  exp  (jn/4 ).  Shown  in 
Figure  P-1  is  the  magnitude  of  this  response  for  a  limited  amplitude  range  of  1 10  dB.  The 
amplitude  range  is  intentionally  limited  to  emphasize  the  shape  of  the  response  both  near  the 
pole  and  near  z  =  0.  It  is  also  useful  to  limit  the  amplitude  range  given  that  the  magnitude 
response  at  the  pole  is  infinite  in  value.  Referring  to  the  figure,  we  note  that  the  response  in  the 
vicinity  of  the  pole  is  very  distinct  and  non-zero.  The  response  evaluated  on  the  unit  circle  at 
z  =  exp  (joS)  is  the  (discrete-time)  Fourier  transfonn,  and  is  of  particular  interest.  Our  given 
pole  response  is  that  of  a  unilateral  complex  sinusoid,  thus  the  peak  on  the  unit  circle  at  the  pole 
is  due  to  the  sinusoidal  nature  of  the  signal  sequence.  The  response  on  the  unit  circle  near  the 
peak  is  partly  due  to  the  sharp  transition  of  the  sequence  as  the  sinusoid  starts.  Various 
definitions  of  signal  bandwidth  exist  to  help  identify  the  frequency  region  of  significant  spectral 
content  about  the  center  frequency  at  the  pole  peak.  We  will  be  able  to  leverage  these 
observations  and  definitions  to  define  a  corresponding  response  of  reduced  bandwidth  after  first 
considering  the  individual  real  and  imaginary  responses  of  our  example.  Shown  in  Figure  P-2  is 
the  magnitude  of  the  imaginary  part  of  the  example  response,  z/(z  —  a).  Once  again,  the 
amplitude  range  is  arbitrarily  chosen  for  plotting  purposes.  It  is  immediately  apparent  that  the 
imaginary  part  of  the  response  contributes  to  the  Fourier  transform  in  the  vicinity  of  the  peak, 
and  therefore  to  the  signal  bandwidth.  This  motivates  us  to  further  consider  the  magnitude  of  the 
real  part  of  the  response  as  presented  in  Figure  P-3.  This  figure  is  particularly  revealing  in  that 
the  unit  circle  is  seen  to  have  an  apparent  constant  low-level  response,  with  the  exception  of  the 
sharp  peak  at  the  pole.  Observations  of  the  actual  numerical  values  of  this  response  on  the  unit 
circle  reveals  that  the  peak  occurs  in  a  very  tight  frequency  region.  Small  (machine  epsilon) 
deviations  from  the  pole  frequency  change  the  response  from  infinity  (not-a-number  or  NaN)  at 
the  pole,  to  a  magnitude  of  0.5  elsewhere.  Motivated  by  this  and  previous  observations 
regarding  the  imaginary  part  of  the  response,  we  consider  defining  a  reduced  bandwidth 
representation  frequency  response  as 

Sa{e^)  =  Re{z/ (z  -  a)} \z=e,a  -  0.5  .  (2-1) 

The  magnitude  of  such  a  response  is  shown  in  Figure  P-4.  Note  that  the  magnitude 
response  on  the  unit  circle  is  distinctly  apparent  and  small  in  value,  and  is  confirmed  to  be  zero 
at  all  frequencies  except  the  exact  value  (within  machine  epsilon  deviations)  of  the  pole 
frequency.  It  is  also  interesting  to  note  that  similar  functions  have  been  sought  for  related  but 
distinct  purposes,  in  defining  non-zero  regions  for  time-frequency  distributions.  The  visual 
appearance  of  these  similar  functions  have  led  authors  to  refer  to  them  as  “bowtie”  functions  for 
obvious  reasons.  (See  e.g.,  [Boashash],  [Brown]).  As  a  distinction,  the  “bowtie”  response 
shown  in  the  figure  is  noted  to  have  an  infinite  peak  at  the  pole  (as  opposed  to  some  finite  value 
sometimes  cited  in  the  references). 
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Figure  P-1.  Magnitude  plot  of  a  single  pole  at  u>  =  rt/4  on  the  unit  circle  of  the  z-plane. 
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Figure  P-2.  Magnitude  plot  of  the  imaginary  component  of  a  single  pole  at  u>  =  tt/4  on  the 
unit  circle  of  the  z-plane. 
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Figure  P-3.  Magnitude  plot  of  the  real  component  of  a  single  pole  at  oo  =  Tt/4  on  the  unit  circle 
of  the  z-plane. 
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Figure  P-4.  Magnitude  plot  of  the  real  component  of  a  single  pole  at  oo  =  tt/4  on  the  unit  circle 
of  the  z-plane,  with  .5  bias  removed. 
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2.1  Mated  Poles 


Having  defined  the  reduced  bandwidth  representation  for  a  single  pole  as  in  Eq.  (2-1),  we 
need  to  identify  the  mathematical/physical  interpretation  of  such  a  representation.  We  begin  by 
defining 

Ya(eJ“)  =  ^(Vw)  +  \Y2{en  ,  (2-2) 


where 


and 


Fi(e'w) 


(2-3) 


Y2(eJt0)  =  Y*(ej0J) 


(2-4) 


is  the  conjugate  of  Vj.  From  Eqs.  (2-3)  and  (2-4),  Eq.  (2-2)  is  equivalent  to  the  first  tenn  of  Eq. 
(2-1).  Substituting  from  Eqs.  (2-3)  and  (2-4)  into  Eq.  (2-2)  and  simplifying  results  in 


Ya(.z)\z=eJ<o  =  0.5 


—  +  1  + 


Replacing  the  first  tenn  of  Eq.  (2-1)  with  the  results  from  Eq.  (2-5)  we  obtain 


(2-5) 


-0.5  .  (2-6) 

z=eia> 

From  Eq.  (2-6)  we  see  that  the  reduced  bandwidth  representation  of  a  single  pole  can  be 
interpreted  as  the  response  of  a  two  pole  system.  The  poles  of  this  new  system  are  identified  as 
the  original  pole,  a,  and  the  conjugate  of  its  reciprocal,  1/a*.  Herein,  these  poles  will  be 
referred  to  as  mated  poles. 

Shown  in  Figure  P-5  is  the  imaginary  component  of  the  mated  pole  response  Ya(z ),  for  a 
specific  example,  a  —  ea+ja)  ,  with  a  —  .  1  and  a>  =  n/\.  Note  that  the  mated  poles  occur  at  the 
same  frequency,  but  at  radial  distances  from  the  origin  that  are  reciprocal  in  magnitude.  This  is 
easily  confirmed  by  the  fact  that  —  =  e~a+]U> .  Figure  P-6  presents  the  magnitude  of  the  Fourier 

transform  of  a  single  pole,  and  that  of  the  associated  reduced  bandwidth  representation.  The 
blue-colored  response  represents  the  original  single  pole  frequency  response,  and  the  red-colored 
response  represents  that  of  the  reduced  bandwidth  response.  For  this  example,  a  =  —0.3  and 
a)  —  7r/4,  with  1024  uniformly  spaced  samples  selected  on  the  unit  circle.  As  seen  in  the  figure, 
the  mated  pole  response  with  bias  removal  is  significantly  reduced  in  bandwidth. 
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Figure  P-5.  Imaginary  component  of  z-plane  response  of  mated  poles  for  a  =  .1  and  oo  =  tt/4 
radians. 
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Figure  P-6.  Fourier  transform  and  Reduced  Bandwidth  Representation  Fourier  transform  for  a 
single  pole  at  u>  =  tt/4  and  a  =  —0.3,  sampled  at  N  =  1024  points  on  the  unit  circle. 
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Another  specific  example  is  presented  in  Figure  P-7.  In  this  case,  a  —  le~3  and 
co  —  n/4,  with  1024  uniformly  spaced  samples  selected  on  the  unit  circle.  Note  that  as  the 
original  pole  approaches  the  unit  circle  both  responses  sharpen  as  expected,  however,  the 
reduction  in  bandwidth  becomes  more  pronounced.  To  gain  further  insights,  we  are  interested  in 
the  inverse  z-transform  of  our  mated  pole  response  identified  by  Eq.  (2-5). 

Consider  the  case  where  a  =  ea+ J0)  is  inside  the  unit  circle,  i.e.,  a  <  0.  The  region  of 
convergence  (ROC)  associated  with  the  term  z/(z  —  a)  for  this  case  is  |z|  >  |a|.  The  ROC  for 
this  term  extends  from  (but  not  including)  a  ring  of  radius  equal  to  the  pole  magnitude,  outwards 
in  all  directions  on  the  z-plane,  and  includes  the  unit  circle.  The  inverse  z-transform  for  this  tenn 
of  the  response  is  well  known,  and  is  equal  to  anit[n].  The  ROC  associated  with  the  term 
— z/(z  —  1/a*)  for  this  case  is  |z|  <  |l/a*|.  It  extends  from  (but  not  including)  a  ring  of  radius 
equal  to  the  magnitude  of  the  pole  1/a*,  inwards  on  the  z-plane  towards  zero,  and  also  includes 
the  unit  circle.  The  inverse  z-transform  for  this  tenn  of  the  response  is  also  well  known,  and  is 
equal  to  anu[—n  —  1].  The  remaining  term  is  a  constant  with  a  ROC  that  includes  the  entire  z- 
plane,  and  is  known  to  have  a  Kronecker  delta  function,  5[n],  as  its  inverse.  Combining  terms, 
the  superposition  principle  allows  us  to  write 

y[n]  —  0.5(anit[n]  +  <5[n]  +  (l/a*)nit[— n  —  1])  .  (2-7) 

Thus  our  mated  poles  result  in  the  sum  of  a  right-sided  sequence  due  to  the  pole  inside  the  unit 
circle,  a  Kronecker  delta  sequence,  and  a  left-sided  sequence  due  to  the  pole  outside  the  unit 
circle.  From  Eq.  (2-7)  and  the  relationship  between  T(e-,w)  and  S(eja> )  given  in  Eqs.  (2-5)  and 
(2-6)  allows  us  to  write 


s[n]  =  0.5(anu[n]  +  (l/a*)nit[— n  —  1])  .  (2-8) 

This  is  our  sequence  resulting  from  the  reduced  bandwidth  representation  of  the  original  signal. 
Such  a  sequence  is  shown  in  Figure  P-8,  with  the  scale  factor  of  Vi  omitted.  The  interpretation  of 
the  result  in  the  sequence  domain  is  that  the  reduced  bandwidth  representation  leads  to  a 
corresponding  sequence  that  exists  for  all  indices,  n.  For  the  example  shown  in  the  figure,  the 
red-  and  magenta-colored  dots  are  respectively  the  real  and  imaginary  parts  of  the  right-sided 
component,  and  the  blue-  and  cyan-colored  dots  are  respectively  the  real  and  imaginary  parts  of 
the  left-sided  component.  In  this  example,  a  =  —.05  and  co  —  n/34.  The  case  where  the  pole, 
a,  is  outside  the  unit  circle  (i.e.  the  original  sequence  is  left-sided)  can  be  analyzed  in  a  similar 
fashion,  but  the  basic  results  are  the  same. 

In  the  sequence  domain,  the  signal  sequence  no  longer  has  an  abrupt  start  near  index  0, 
thereby  reducing  potentially  extraneous  bandwidth  in  the  Fourier  domain.  Also,  the  new 
sequence  is  seen  to  be  symmetric  in  magnitude. 
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Figure  P-7.  Fourier  transform  and  reduced  bandwidth  representation  Fourier  transform  for  a 
single  pole  at  u>  =  tt/4  and  a  =  le-3,  sampled  at  N  =  1024  points  on  the  unit  circle. 


Example  Sequence  from  a  Reduced  Bandwidth  Representation 
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Figure  P-8.  An  example  sequence  resulting  from  a  reduced  bandwidth  representation  for  a 
decaying  (right-sided  sequence)  with  a  =  —.05  at  a  frequency  of  oo  =  rr/34  radians. 
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3.0  Application  and  Results 


Before  presenting  results  from  various  test  signal  inputs,  a  generalization  of  the  reduced 
bandwidth  representation  is  required.  The  Fourier- transform  of  our  input  sequence  segment, 
x[n],  can  be  decomposed  as 

X{eJ”)  =  Xd(e&)  +  Xg(en  ,  (3-1) 

where 

Xd(e^)=YliAiXi(en  (3-2) 

is  the  sum  of  the  Fourier-transforms  of  each  right-sided  component,  Xj[n]  =  cqnit[n],  and 

Xg(e*)=XkAkXk(en  (3-3) 

is  the  sum  of  the  Fourier-transforms  of  each  left-sided  component,  xk[n ]  =  aknu[—n  —  1],  The 
(complex  valued)  factors,  AL  and  Ak,  are  the  amplitudes  of  the  right-sided  and  left-sided 
components  respectively.  From  Eq.  (2-5)  along  with  Eqs.  (3-1)  through  (3-3),  it  can  be  shown 
that  the  Fourier  transform  resulting  from  the  combined  mated  pole  responses  is 

y{e>“)  =  E(M»eWe/")})  -  JWKepWe'")})  •  (3-4) 

Finally,  corresponding  to  Eq.  (2-1)  we  obtain  the  reduced  bandwidth  Fourier  representation 

S(eJco)  =  XiAiCRelz/Cz  -  af)} \z=e]a>  -  0.5)  -  '£kAk(Re{z/(z  -  ak)} \z=eJo>  -  0.5)  .  (3-5) 

It  should  be  noted  that  in  deriving  Eq.  (3-5),  we  have  ignored  any  sequence  shift  associated  with 
each  of  the  left-sided  components  that  arises  due  to  the  length  of  the  segment,  L,  being  analyzed. 
Unless  otherwise  specified,  application  results  presented  herein  ignore  the  effects  of  this  shift. 

For  some  uses  of  the  reduced  bandwidth  Fourier  representation,  taking  this  shift  into  account 
may  be  appropriate.  This  may  be  particularly  true  when  considering  implementations  of  filtering 
in  the  Fourier  domain  via  products  of  transforms.  Detennining  the  potential  for  such  use  is 
reserved  for  future  research.  For  this  report,  the  factor  required  to  account  for  this  shift  is 
omitted  to  avoid  ripple  artifacts  in  magnitude  spectra  that  can  be  visually  distracting. 

3.1  Signal  Segment  Decomposition 

To  make  use  of  the  above  results,  a  method  of  decomposing  an  input  signal  data  segment 
into  component  right-  and  left-sided  exponential  sequences  is  required.  Relevant  methods  do 
exist  and  are  found  in  the  literature.  However,  the  methods  typically  decompose  into  decaying 
right-sided  exponentials  only,  or  decompose  into  decaying  and/or  growing  right-sided 
exponentials.  One  such  method  is  the  Matrix  Pencil  (MP)  algorithm  (See  e.g.  [Sarkar]),  and  is 
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adopted  for  use  in  this  research.  The  MP  algorithm  is  attractive  in  that  it  directly  yields  most  of 
what  is  required  to  estimate  the  reduced  bandwidth  Fourier  representation.  In  particular,  the  MP 
algorithm  will  generate  estimates  of  the  parameters  { |  d  t| ,  (pL,  ,  eo,}  for  the  decaying  and/or 
growing  right-sided  complex  exponentials.  These  parameters  are  the  component  magnitude, 
phase  (in  radians),  decay  (or  growth)  and  radian  frequency  respectively.  Note  that  the  complex 
valued  amplitude  can  be  also  represented  as  At  =  For  the  growing  right-sided  complex 

exponentials,  i.e.,  the  components  where  aL  >  0,  the  conversion  required  is 

Ak  =  C Ai/ade a*L  =  Aiea^  .  (3-6) 

Eq.  (3-6)  is  applied  to  the  amplitude  parameters  generated  by  the  MP  algorithm,  for  the  growing 
components.  The  remaining  parameters  are  unchanged,  and  are  associated  with  the 
corresponding  converted  component  as  indexed  by  the  variable  k. 

The  MP  algorithm  allows  for  user-selection  of  the  total  number  of  components,  Kc,  and 
allows  for  up  to  L/ 2  components  to  be  specified.  For  our  purposes,  a  variant  of  the  MP 

algorithm  has  been  created  that  allows  for  user  specification  of  up  to  (j)  —  1  components,  to 

help  mitigate  instabilities  that  can  otherwise  occur.  These  instabilities  are  detected  and  then 
avoided  by  increasing  the  specified  number  of  components  by  1,  as  used  in  [Haddad]. 

Thus  we  can  use  the  Matrix  Pencil  algorithm  for  signal  segment  decomposition  and  use 
the  associated  parameter  estimates  to  generate  a  reduced  bandwidth  Fourier  representation  for 
consecutive  signal  segments.  The  resulting  time-frequency  representation  is  referred  to  herein  as 
the  pencil  gram. 

3.1.1  Decomposition  Iterations 

To  mitigate  concerns  regarding  the  accuracy  of  signal  models  and  decompositions 
employed  in  this  research,  a  simple  iteration  technique  has  been  developed.  Unless  otherwise 
specified,  this  iteration  technique  is  applied  to  the  multi-segment  examples  presented  herein.  In 
particular,  a  two-pass  iteration  is  used  to  address  model  inaccuracies  as  a  result  of  limitations  on 
the  number  of  signal  components,  Kc. 
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3.2  Example  Pencilgrams 

While  specific  metrics  could  be  introduced  to  assess  the  performance  of  any  given  time- 
frequency  representation,  such  metrics  tend  to  be  relevant  only  to  a  subset  of  signal  scenarios. 
For  example,  one  could  consider  single-tone  input  tests  for  assessing  frequency  resolution  and 
second-order  hannonic  distortion.  Likewise,  two-tone  tests  could  be  devised  to  assess  third- 
order  distortions.  Such  objective  tests  of  the  pencilgram  and  other  reduced  bandwidth  Fourier 
representations  are  reserved  for  future  research.  In  lieu  of  specific  metrics,  the  signal  processing 
community  has  often  relied  on  subjective  analyses  of  signals  of  interest  when  evaluating  time- 
frequency  representations.  To  demonstrate  the  potential  benefits  of  the  pencilgram,  the  latter 
approach  is  taken  for  this  report. 

3.2.1  Kay  and  Marple  Test  Signal  Segment 

First,  we  can  leverage  published  subjective  comparisons  as  applied  to  a  specially 
constructed  test  signal  proposed  in  [Kay].  This  signal  consists  of  a  band-limited  (band-pass) 
noise  process  centered  at  0.7n  radians  plus  three  tones,  each  at  0.27T,  0.47T  and  0.427T  radians. 
The  amplitudes  of  the  sinusoids  were  0.1,  1.0  and  1.0  giving  rise  to  +10,  +30  and  +30  dB  signal- 
to-noise  power  ratios  respectively.  The  segment  length  is  L  —  64.  (No  other  information  was 
given  regarding  the  filter  used  to  generate  the  band-pass  noise.)  Shown  in  Figure  KM-1  is  the 
Fast  Fourier-Transform-  (FFT-)  based  spectral  estimate  of  the  test  signal.  The  samples  were  not 
weighted  prior  to  the  FFT  (i.e.,  a  rectangular  window  was  used),  nor  was  zero-padding 
employed.  Because  of  the  short  segment  length,  frequency  resolution  is  limited  as  seen  in  the 
figure.  In  particular,  although  the  first  tone  is  distinguishable  as  a  component,  the  two  remaining 
tones  are  too  close  in  proximity  to  be  able  to  distinguish  them  as  separate  signal  components. 

The  energy  due  to  the  band-pass  noise  is  also  apparent,  although  it  would  be  difficult  from  this 
figure  to  identify  this  component  as  band-pass  noise;  it  could  erroneously  be  identified  as  an 
additional  four  or  five  sinusoids.  This  may  be  related  to  the  order  and  type  of  filter  used  to  filter 
the  noise.  For  example,  a  high  order  (relative  to  L )  finite  impulse  response  (FIR)  filter  would 
require  a  longer  segment  to  allow  the  FFT  to  give  a  more  accurate  representation  of  the  filter 
shape.  This  fact  combined  with  the  nature  of  the  noise  signal  itself,  would  not  lead  one  to  expect 
to  see  the  filter  shape  alone  as  implied  in  the  reference,  but  rather  the  product  of  the  filter 
response  and  the  true  underlying  noise  spectrum  for  a  given  segment  length.  Thus  at  best,  we 
would  only  expect  to  see  something  similar  to  what  is  observable  in  the  figure  for  the  band-pass 
noise.  Note  in  particular  that  the  FFT-based  spectral  estimate  does  fairly  well  with  regard  to  the 
strengths  of  the  sinusoids;  they  are  commensurate  with  the  known  signal  strengths. 
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For  comparison,  the  (single  segment)  pencilgram  of  the  same  sequence  is  shown  in 
Figure  KM-2.  For  this  pencilgram,  18  poles  are  requested  and  the  frequency  response  is 
generated  via  5000  unifonnly  spaced  samples  on  the  unit  circle.  Sharp  peaks  are  noted  at  the 
locations  of  the  true  sinusoids,  and  the  shape  of  the  band-pass  noise  is  apparent.  Note  that  the 
magnitude  range  on  this  figure  is  50  dB  greater  than  that  of  the  previous  figure.  Overall,  the 
pencilgram  for  this  sequence  compares  favorably  with  the  FFT-based  method,  however,  it  is 
noted  that  the  magnitudes  of  the  sinusoidal  components  are  much  larger  than  that  of  the  FFT- 
based  method.  It  is  also  interesting  to  note  that  increasing  the  number  of  requested  poles  to  3 1 
further  sharpened  the  response  due  to  the  sinusoids,  and  also  increased  their  associated 
magnitudes.  In  general,  as  poles  approach  the  unit  circle,  they  give  rise  to  large  amplitudes  in 
very  small  frequency  ranges. 
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Figure  KM-1.  An  FFT-based  spectral  estimate  of  the  test  signal  identified  in  [Kay];  rectangular 
data  window  used  and  no  zero  padding. 
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Figure  KM-2.  Pencilgram  of  the  test  signal  identified  in  [Kay];  Kc  =  18  complex  exponential 
components  requested  with  5000  samples  on  the  unit  circle. 

3.2.2  Sinusoid  and  Growing  Right-Sided  Exponential  Signal  Segment 

The  next  example  consists  of  a  unit  amplitude  sinusoid  at .  2n  radians  plus  a  single 
complex  exponential  component.  The  complex  exponential  has  an  amplitude  of  1  (prior  to 
amplitude  conversion),  with  a  =  .09  at  a  frequency  of  oj  =  1.4  radians.  Note  that  the  left-sided 
representation  of  this  signal  has  a  frequency  of -1.4  radians,  which  will  peak  at  2n  —  1.4  radians, 
due  to  the  modulo  nature  of  the  frequency  response.  The  FFT-based  spectral  estimate  is 
presented  in  Figure  KM-3.  From  this  figure  it  is  apparent  that  because  of  the  overwhelming 
strength  of  the  exponential  component  relative  to  the  sinusoid,  the  sinusoid  itself  is  barely 
detectable  as  a  small  glitch  in  the  spectrum  at  the  expected  locations  of .  2n  and  2n  —  .2n  — 

1.8n  radians. 

The  (single  segment)  pencilgram  of  this  same  sequence  is  presented  in  Figure  KM-4. 

Both  components  are  clearly  distinguishable.  We  once  again  note  that  the  pole  estimate  that  is 
near  the  unit  circle  due  to  the  sinusoidal  component,  leads  to  large  amplitude  peaks 
approximating  Dirac  delta  functions.  Thus,  the  function  identified  in  Eq.  (2-1)  can  be  considered 
a  nascent  (Dirac)  delta  function.  More  will  be  discussed  later  in  this  report  regarding  this 
property. 
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3.2.3  Multi-Segment  Sinusoid 

We  now  progress  to  test  signals  that  are  more  substantial  in  length,  and  analyzed  by 
partitioning  the  sequence  into  consecutive  overlapping  segments  of  short  length.  The  overlap  for 
all  subsequent  examples  is  50%,  and  segment  lengths  are  L  —  64  samples.  It  should  also  be 
noted  that  these  pencilgrams  have  been  generated  as  a  result  of  a  2-pass  Matrix  Pencil  method 
described  previously.  The  number  of  unifonnly  spaced  samples  on  the  unit  circle,  N,  is  specified 
as  1024  for  the  pencilgrams.  All  spectrograms  presented  for  comparisons  were  calculated  using 
Hanning  window  weighting  and  were  zero-padded  to  length  N=  1 024. 

The  first  such  example  is  a  pure  sinusoid  at  u>  =  66n/N  radians,  and  is  shown  in  Figure 
S-l.  The  upper  subplot  is  the  pencilgram  and  the  lower  subplot  is  the  FFT-based  spectrogram. 
For  comparison  purposes,  the  amplitude  (magnitude)  range  on  both  subplots  was  set  from  -40  to 
+40  dB.  Careful  inspection  of  the  pencilgram  shows  a  sharp  line  at  the  correct  frequency.  All 
signal  energy  is  contained  at  this  frequency  index,  and  all  others  contain  no  measurable  energy. 
In  contrast,  the  spectrogram  has  substantial  energy  spread  in  FFT  frequency  bins  at  and  around 
the  correct  frequency.  Also,  the  spectrogram  contains  signal  energy  in  all  other  frequency  bins 
resulting  from  assumed  segment  periodicity  and  the  implicit  convolution  of  the  frequency 
response  of  the  window  with  that  of  the  data.  Figure  S-2  shows  the  spectrogram  of  the  original 
signal  and  the  spectrogram  of  the  error  associated  with  the  MP  signal  model.  As  seen  in  figure, 
the  signal  model  for  this  example  is  highly  accurate. 
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Figure  KM-3.  An  FFT-based  spectral  estimate  of  a  test  signal  consisting  of  a  unit  amplitude 
sinusoid  at  a)  =  ,2n  radians  plus  a  unit  amplitude  growing  (right-sided)  complex  exponential 
with  a  =  .09  and  frequency  of  1.4  radians;  rectangular  data  window  used  and  no  zero  padding. 
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Figure  KM-4.  Pencilgram  of  a  test  signal  consisting  of  a  unit  amplitude  sinusoid  at  a)  =  .2tt 
radians  plus  a  unit  amplitude  growing  (right-sided)  complex  exponential  with  a  =  .09  and 
frequency  of  1.4  radians;  Kc  =  18  complex  exponential  components  requested  with  N  =  5000 
samples  on  the  unit  circle. 
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2-pass  //  PEncilGram  (tap)  S  Spectrogram  (bottom);  File:  Sinusoid 


Tima  Segment 


Figure  S-l.  Pencilgram  (upper  subplot)  and  Spectrogram  (lower  subplot)  of  a  sinusoidal  input 
test  signal;  2x28=56  complex  exponential  components  requested  for  the  pencilgram;  frequency 
bins  0  to  100  shown  and  color-map  chosen  for  visibility. 


2-pass  //  Spectrogram  (D  to  .5N-I)  B  Error  Spectrogram  (.5N  to  N-l);  File:  Sinusoid 
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Figure  S-2.  Spectrogram  of  error  signal  (bins  ,5N  to  N-l)  and  spectrogram  of  sinusoidal  input 
test  signal  (bins  0  to  .5N-1);  2x28=56  complex  exponential  components  requested  for  the 
pencilgram. 
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3.2.4  Multi-Segment  Chirp 

For  this  example,  a  chirp  signal  (a  single  swept  frequency)  is  used  to  generate  both  the 
pencilgram  and  FFT-based  spectrogram  as  shown  in  the  corresponding  subplots  of  Figure  C-l. 
Once  again,  the  amplitude  (magnitude)  range  is  set  from  -40  to  40  dB  for  comparisons.  The 
modulation  index  is  intentionally  chosen  to  be  large  to  concentrate  energy  near  the  instantaneous 
frequency.  Referring  to  the  spectrogram,  while  we  expect  some  energy  spread  near  the  peak 
energy,  again  the  implicit  segment  periodicity  and  data  window  convolution  produces  significant 
spectral  artifacts  throughout  the  spectrum.  In  contrast,  the  pencilgram  shows  significant  energy 
near  the  peak  and  quickly  subsides  outside  this  location.  This  is  true  for  all  time  segments. 
Looking  across  time  segments,  a  variation  is  noted  in  this  spread  of  energy  away  from  the  peak. 
However,  it  is  also  pointed  out  that  the  level  of  this  energy  (shown  in  bright  blue)  is  about  30  dB 
lower  than  that  of  the  spectrogram.  This  is  an  example  where  a  more  detailed  analysis  could  be 
performed  to  obtain  objective  metrics  regarding  this  apparent  spectral  spread.  In  particular,  the 
true  spectrum  determined  analytically  could  be  compared  to  the  results  obtained  in  the 
pencilgram.  Such  objective  analyses,  now  motivated  by  the  subjective  results,  are  reserved  for 
future  research.  Shown  in  Figure  C-2  is  the  spectrogram  of  both  the  error  due  to  MP  modeling 
inaccuracy  and  the  original  signal.  Again  note  that  the  MP  model  is  quite  accurate  as  indicated 
by  the  level  of  the  error  signal  relative  to  the  original.  This  is  accomplished  in  this  example  with 
14  components  requested  per  iteration,  leading  to  potentially  2(1  +  14)  =  30  complex  sinusoids 
generated  by  the  MP  algorithm.  Note  that  up  to  32  components  could  have  been  generated  by 
the  MP  algorithm  for  each  64  sample  segment  per  iteration.  Thus  for  two  iterations  as  used  in 
this  example,  64  complex  sinusoids  and  associated  parameters  could  have  been  generated. 

As  previously  described,  the  phase  factor  e~](x)L  has  been  omitted  in  our  development  of 
the  reduced  bandwidth  representation.  Comparing  the  results  shown  in  Figures  C-3  and  C-4,  we 
observe  that  incorporation  of  this  factor  changes  the  corresponding  pencilgram.  Figure  C-3 
shows  the  pencilgram  without  this  factor,  and  Figure  C-4  shows  the  pencilgram  incorporating 
this  factor  for  all  left-sided  components  in  the  decomposition.  The  merit  in  incorporating  this 
phase  factor  is  left  for  further  study;  this  analysis  is  expected  to  benefit  from  the  development  of 
the  objective  metrics  already  alluded  to. 

Results  from  a  variation  of  the  chirp  test  signal  are  also  presented  here.  This  variation  is 
a  band-pass  filtered  version  of  the  previously  described  chirp  signal.  The  band-pass  filter  used  to 
generate  this  signal  is  an  Butterworth  infinite  impulse  response  (HR)  filter  as  given  by  the 
Matlab  command  [bfilt,  afilt]  =butter  (16,  [ .  4  .  6] )  .  This  command  generates  a  16th 
order  band-pass  filter  centered  at  the  middle  of  the  Nyquist  band,  with  200  dB  attenuation  at .  2n 
and  .  87T  radians,  and  over  400  dB  attenuation  at  the  band  edges.  Figures  FC-1  and  FC-2  show 
the  results.  Note  that  the  pencilgram  gives  intuitively  satisfying  results,  whereas  the  spectrogram 
again  has  significant  artifacts  that  obfuscate  the  true  spectral  content. 
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2-pass  //  PencilGram  (top)  &  Spectrogram  (bottom);  File:  Chirp 
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Figure  C-l.  Pencilgram  (upper  subplot)  and  Spectrogram  (lower  subplot)  of  a  chirp  (swept 
frequency)  input  test  signal;  2x28=56  complex  exponential  components  requested  for  the 
pencilgram. 
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Figure  C-2.  Spectrogram  of  error  signal  (bins  .5N  to  N-1)  and  spectrogram  of  chirp  (swept 
frequency)  input  test  signal  (bins  0  to  .5N-1);  2x28=56  complex  exponential  components 
requested  for  the  pencilgram. 
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2-pass  //  PencilGram;  File:  Chirp 


Time  Segment 


Figure  C-3.  Pencilgram  of  a  chirp  (swept  frequency)  input  test  signal;  2x28=56  complex 
exponential  components  requested  for  the  pencilgram. 


2-pass  //  PencilGram;  File:  Chirp 
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Figure  C-4.  Pencilgram  of  a  chirp  (swept  frequency)  input  test  signal;  2x28=56  complex 
exponential  components  requested  for  the  pencilgram.  Left  sided:  X  =  X  -  exp(- 
j*W*L)*exp(alphak*(L-l))*exp(j*Wk*(L-l))*r(k)*exp(j*ph(k)).*F  (See  Appendix). 
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2-pass  //  PencilGram  (top)  &  Spectrogram  (bottom);  File:  Band-pass  Filtered  Chirp 
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Figure  FC-1.  Pencilgram  (upper  subplot)  and  Spectrogram  (lower  subplot)  of  a  band-pass 
fdtered  chirp  (swept  frequency)  input  test  signal;  2x28=56  complex  exponential  components 
requested  for  the  pencilgram. 
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Figure  FC-2.  Spectrogram  of  error  signal  (bins  ,5N  to  N-l)  and  spectrogram  of  filtered  chirp 
(swept  frequency)  input  test  signal  (bins  0  to  .5N-1);  2x28=56  complex  exponential  components 
requested  for  the  pencilgram. 
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3.2.5  Multi-Segment  Speech  Signal 

The  last  results  presented  in  this  section  are  obtained  for  a  specific  speech  signal  selected 
from  the  TIMIT  database,  file  faksO_sal.wav.  This  is  a  fairly  clean  speech  signal  of  a  short 
utterance  (3  or  4  seconds)  from  a  female  speaker,  sampled  at  8000  samples  per  second.  In 
general,  speech  represents  a  “spectrally  diverse”  signal  containing  many  components  including 
some  that  are  both  narrow-  and  wide-band  nature.  As  a  spectrally  rich  signal,  it  also  serves  to 
demonstrate  the  utility  of  the  2-pass  iteration  technique  previously  developed  and  used  in  this 
research  to  enhance  MP  model  and  signal  decomposition  accuracy. 

In  Figures  SP-la  and  SP-lb,  we  demonstrate  signal  model  and  decomposition 
enhancement  obtained  via  2-pass  iteration.  In  the  first  figure,  we  show  the  spectrograms  of  both 
the  original  speech  (lower  portion)  and  the  model  error  (upper  portion)  when  two  iterations  of 
decomposition  are  performed.  It  is  noted  that  relative  to  using  one  iteration  as  seen  in  Figure  SP- 
lb,  the  model  error  is  significantly  reduced  by  the  second  iteration. 

For  comparative  analyses,  the  pencilgram  and  spectrogram  of  the  selected  speech  signal 
are  presented  in  Figure  SP-2.  The  amplitude  range  for  these  plots  have  been  set  from  -70  to  +10 
dB.  Observe  that  the  pencilgram  provides  significantly  enhanced  frequency  resolution  relative  to 
the  FFT-based  spectrogram.  This  is  particularly  helpful  in  the  analyses  of  speech  signals,  given 
the  harmonic  content  present  in  the  signal.  Speech  waveforms  can  be  described  as  having 
voiced,  unvoiced  and  silence  regions,  along  with  transitions  between  these  regions,  as  the  signal 
progresses  in  time.  The  voiced  regions  and  corresponding  segments  tend  to  have  this  harmonic 
structure.  These  regions  are  of  particular  interest  to  researchers  for  a  variety  of  applications 
including  automated  speaker  identification  and  automated  speech  recognition. 

In  general,  the  amplitude  range  for  any  time-frequency  representation  must  be  carefully 
selected,  and  the  more  useful  ranges  are  test  signal  dependent.  A  new  amplitude  range  of  -40  to 
+40  dB  has  been  selected  for  the  pencilgram  presented  in  Figure  SP-3.  It  can  be  observed  in  this 
and  the  previous  figure,  the  pencilgram  exhibits  spectral  energy  across  a  large  portion  of  the 
spectrum  for  some  time  segments  (e.g.,  between  370  and  400).  Conversely,  in  the  FFT-based 
spectrogram,  these  wide-band  events  do  not  appear  as  prominently.  The  question  then  arises  as 
to  the  significance  of  these  events,  and  whether  or  not  they  accurately  portray  the  true  underlying 
spectral  content. 

One  possibility  is  that  the  potential  inaccuracy  is  a  result  of  the  omitted  phase  factor, 
e]a)L,  for  each  left-sided  component.  Initial  comments  can  be  made  regarding  this  possibility  by 
referring  to  Figure  SP-4.  In  this  figure,  the  pencilgram  has  been  generated  with  the  inclusion  of 
this  factor.  Referring  to  the  figure,  although  there  are  amplitude  ripples  in  the  frequency 
dimension,  the  wide-band  events  are  still  present. 
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Figure  SP-la.  2-pass  Matrix  Pencil  decomposition;  spectrogram  of  error  signal  (bins  .5N  to  N- 
1)  and  spectrogram  of  a  speech  input  test  signal  (bins  0  to  .5N-1);  Kc  =  28  complex  exponential 
components  requested  for  the  decomposition  for  each  iteration. 
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Figure  SP-lb.  1-pass  Matrix  Pencil  decomposition;  spectrogram  of  error  signal  (bins  .5N  to  N- 
1)  and  spectrogram  of  a  speech  input  test  signal  (bins  0  to  .5N-1);  Kc  =  28  complex  exponential 
components  requested  for  the  decomposition. 
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Figure  SP-2.  Pencilgram  (upper  subplot)  and  spectrogram  (lower  subplot)  of  a  speech  input  test  signal; 
2x28=56  complex  exponential  components  requested  for  the  pencilgram. 
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Figure  SP-3.  Pencilgram  of  a  speech  input  test  signal;  2x28=56  complex  exponential 
components  requested  for  the  pencilgram. 
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Regardless  of  the  accuracy  issue,  it  is  interesting  to  note  that  these  events  can  be 
selectively  diminished  or  removed  in  the  pencilgram  as  a  direct  characteristic  of  the  MP 
algorithm.  Because  each  component  resulting  from  the  MP  algorithm  is  parameterized  as 
previously  described,  those  components  having  large  values  of  |  a  |  can  be  detected  and  handled 
differently  than  other  components.  For  example,  only  those  components  with  values  of  |cr|  less 
than  some  specified  value  can  be  used  to  generate  the  pencilgram.  This  is  the  case  shown  in 
Figure  SP-5.  Here,  components  with  values  less  than  .08  in  magnitude  were  used.  Note  that  this 
has  a  significant  effect  in  the  wide-band  event  regions  previously  identified  (e.g.,  segments  370 
to  400). 


In  general,  a  wide  variety  of  pencilgrams  could  be  conceived  and  presented.  Rather  than 
component  selection  based  on  a  value,  any  other  parameter  or  combination  of  parameters  could 
be  used  to  cue  other  weightings  or  handling  of  the  associated  components. 
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2-pass  //  PencilGram;  File:  faksO  al.wav 
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Figure  SP-4.  Pencilgram  of  a  speech  input  test  signal;  2x28=56  complex  exponential 
components  requested  for  the  pencilgram.  Left  sided:  X  —  X  -  exp(-j*W*L)*exp(alphak*(L- 
l))*exp(j*Wk*(L-l))*r(k)*exp(j*ph(k)).*F  (see  Appendix). 


2-pass  //  PencilGram;  File:  faksOsa1.wav 
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Figure  SP-5.  Pencilgram  of  a  speech  input  test  signal;  2x28=56  complex  exponential 
components  requested  for  the  pencilgram;  a  magnitudes  less  than  .08  only. 
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4.0  Summary,  Discussions  and  Future  Work 


Motivated  by  the  effectiveness  of  the  Matrix  Pencil  (MP)  algorithm  in  decomposing 
short-term  but  otherwise  spectrally  rich  speech  segments,  a  new,  high-resolution  time-frequency 
representation  has  been  conceived.  Background  on  the  application  of  the  MP  algorithm  to 
speech  is  given  in  [Haddad].  The  MP  algorithm  is  characterized  as  a  method  of  decomposing  an 
input  signal  segment  into  constituent  decaying  and/or  growing  complex  exponentials  in  the 
sequence  domain,  with  corresponding  poles  in  the  z-transfonn  domain.  Knowledge  of  the 
relationship  between  the  z-domain  and  the  Fourier  domain  lends  further  credence  to  the  potential 
for  a  new  Fourier  representation.  Investigations  into  the  nature  of  z-domain  poles  in  general, 
leads  to  a  method  of  achieving  a  highly  desirable  characteristic  in  the  new  representation.  This 
characteristic  is  referred  to  herein  as  a  reduced  bandwidth  representation.  As  introduced  in  this 
report,  by  properly  augmenting  the  signal  model  using  additional  poles,  component  bandwidths 
are  significantly  reduced.  In  the  sequence  domain,  this  augmented  signal  model  is  describable  as 
a  conjugate-symmetric  extension  of  the  samples  within  any  given  signal  segment.  This  effective 
window  expansion  mitigates  adverse  effects  of  data  windowing  by  more  accurately  modeling 
many  practical  signals.  One  important  example  of  where  signal  model  enhancement  occurs  is 
for  segmented  sinusoidal  components.  The  enhanced  model  not  only  allows  for  modeling  the 
components  as  eternal  sinusoids,  it  also  allows  for  a  corresponding  representation  in  the  Fourier 
domain  that  is  essentially  a  Dirac  delta  function  as  predicted  by  theory.  Conversely,  existing 
techniques  either  model  all  components  of  a  segment  as  sinusoidal,  or  do  not  yield  infinitesimal 
bandwidths  in  the  Fourier  domain  for  sinusoidal  components.  Moreover,  without  reduced 
bandwidth  representation,  the  remaining  complex  exponential  components  would  otherwise 
result  in  wide-band  content  in  the  Fourier  domain  that  interferes  with  the  visualization  and 
potential  detection  of  other  signal  components. 

In  addition  to  presenting  the  mathematical  development  of  the  reduced  bandwidth 
representation,  a  specific  application  is  also  described.  The  MP  algorithm  has  been  used  with 
this  new  representation  to  create  what  we  refer  to  as  the  pencilgram,  named  in  recognition  of  the 
algorithm.  As  with  many  time-frequency  tools,  subjective  perfonnance  evaluations  have  been 
performed.  Results  have  been  presented  for  a  variety  of  stressing  waveforms,  including  tones, 
swept  tones  and  a  representative  speech  signal.  The  achieved  frequency  resolution  at 
unprecedentedly  small  segment  lengths  as  shown  in  the  examples  clearly  justifies  both  use  of  and 
further  investigations  into  the  properties  of  the  pencilgram. 

4.1  Existing  Signal  Extension  Models 

Techniques  that  model  the  signal  outside  the  extent  of  the  analysis  window  do  already 
exist  in  the  literature.  The  most  well  known  but  sometimes  overlooked  is  that  of  the  DFT.  The 
DFT  is  describable  as  a  model-based  method,  and  as  a  related  algorithm,  the  same  is  true  for  the 
FFT.  The  DFT  models  the  samples  in  the  signal  segment  as  the  weighted  summation  of  an 
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orthogonal  set  of  sinusoids.  However,  it  also  models  the  overall  signal  as  consisting  of  repeated 
consecutive  segments  that  are  identical  to  the  segment  being  analyzed.  Thus  the  analysis 
window  is  in  effect,  expanded  to  include  all  other  samples.  The  problem  is  the  accuracy  of  such 
an  expansion  model.  It  is  rarely  the  case  that  the  segment  being  analyzed  is  extracted  from  a 
sequence  that  is  periodic  with  period  L  samples  (Z  being  the  number  of  samples  in  the  segment). 

In  [Kay],  a  spectral  estimate  is  discussed  that  models  exponential  decaying  components 
outside  the  analysis  window  and  is  referred  to  in  the  publication  as  the  “symmetric  envelope 
exponential  model”.  The  corresponding  extended  signal  x(t)  is  presented  as  a  continuous-time 
signal  comprised  of p  components,  and  the  corresponding  Fourier  transform  is  given  in  the 
reference  as 


X(f)  =  Zm= 1  An  exp  (JGm) 


_ ^am _ 

“4+(2i[/-y2  ’ 


(4-1) 


with  component  amplitudes  Am,  phases  0m  (in  radians),  damping  factors  am,  and  frequencies  fm 
(in  Hz).  The  parameters  of  the  component  exponentials  are  found  via  the  Prony  method  as 
described  in  the  paper.  The  paper  describes  Eq.  (4-1)  as  being  one  possible  Prony  spectrum  and 
is  generated  under  the  assumption  that  all  damping  factors  are  negative  so  that  the  components 
are  decaying  exponentials. 


This  is  in  contrast  to  our  sampled  signal  representation  which  accounts  for  both  decaying 
and  growing  components  as  previously  described.  For  comparison,  corresponding  to  each  term 
of  Eq.  (4-1),  the  reduced  bandwidth  representation  simplifies  from  Eqs.  (2-1)  and  (3-5)  to 


Sai(eJ“)  =  0.5At 


2  cos(o)— o)j)— (e“i+e  ai) 


(4-2) 


for  right-sided  (decaying  with  increasing  index)  components  and 


2  cos(o)— a»j)— (e“i+e  ai) 


(4-3) 


for  left-sided  (growing  with  increasing  index)  components.  These  terms  would  be  summed  over 
both  i  and  k  to  obtain  the  overall  reduced  bandwidth  Fourier  representation,  S(eJ0)).  Note  that 
the  frequency,  m,  is  in  normalized  radians,  and  can  be  converted  to  Hz  with  knowledge  of  the 
sampling  rate,  Fs,  used  to  acquire  a  corresponding  continuous-time  signal.  Specifically,  the 
conversion  is  =  2nf/Fs  .  (See  also  [Oppenheim].)  It  should  also  be  pointed  out  for  clarity  that 
some  acquired  signals  have  been  translated  in  frequency  either  prior  to  or  as  a  result  of  the 
sampling  process.  This  also  should  be  taken  into  account  for  correct  interpretation  of  converted 
frequencies. 
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4.2  Future  Work 


A  special  issue  that  occurs  with  the  reduced  bandwidth  representation  is  the  need  to 
properly  handle  the  nascent  Dirac  delta  functions  that  can  arise.  Having  such  delta  functions 
creates  the  problem  of  determining  how  to  display  segments  with  such  components  on  the  same 
display  as  all  other  segments.  Another  related  problem  is  how  to  associate  the  weights  of  these 
delta  functions  with  the  displayed  amplitudes.  Various  ways  could  be  conceived  to  handle 
nearly  sinusoidal  components  as  special  cases.  A  rudimentary  way  of  addressing  this  issue  is  to 
use  the  FFT  of  the  segment  to  limit  the  amplitudes  of  the  reduced  bandwidth  representation.  An 
example  of  this  is  presented  in  Figure  F-l.  In  this  example,  the  segment  is  weighted  with  a 
Hanning  window,  and  then  zero-padded  to  N  samples  in  length.  For  comparison,  see  also  Figure 
KM-2.  Note  that  most  of  the  pencilgram  remains  unchanged,  with  the  exception  of  the  tone 
components  located  at  0.27T,  0.47T  and  0.427T  radians.  In  addition  to  loss  in  resolution,  a 
disadvantage  of  this  approach  is  the  potential  for  undesired  reductions  in  amplitude  for  otherwise 
low  level  regions  as  seen  in  the  figures. 
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Figure  F-l.  An  example  modified  pencilgram  of  the  test  signal  identified  in  [Kay];  Kc  =  18  complex 
exponential  components  requested  with  5000  samples  on  the  unit  circle. 


30 


The  subjective  results  presented  in  this  research  serve  as  motivation  to  further  pursue 
such  algorithm  refinements  and  to  define  objective  performance  metrics  for  applications  of 
interest.  Objective  performance  measures  will  generally  be  application  specific.  For  example, 
for  potential  application  in  co-channel  interference  mitigation,  a  specific  goal  may  be  to  detect 
and  remove  undesirable  tones  from  the  underlying  signal  of  interest.  Such  tones  can  have 
adverse  effects  in  demodulation  and  other  subsequent  signal  processing  steps.  In  this  example, 
signal-to-interference  (S/I)  ratios  along  with  detection  and  false  alarm  rates  can  serve  as  metrics 
to  determine  if  representations  such  as  the  pencilgram  can  be  used  to  enhance  detection 
performance.  Based  on  interests  in  methods  of  reassignment  as  in  [Gardner],  the  pencilgram 
may  be  useful  as  the  first  step  in  such  techniques,  rather  than  starting  with  e.g.,  the  FFT.  Other 
potential  application  examples  include  Fourier  domain  filtering,  signal  pre-conditioning, 
automated  speech  and  speaker  identification  and  various  pattern  recognition  problems. 

Likewise,  appropriate  metrics  can  be  defined  to  measure  possible  performance  enhancements  for 
these  applications.  It  should  also  be  noted  that  by  identifying  appropriate  metrics,  the  merit  of 
using  other  methods  of  signal  decomposition  could  be  investigated.  In  particular,  it  is  pointed 
out  that  models  that  allow  for  higher  order  signal  components  could  be  investigated.  Techniques 
to  utilize  such  models  may  be  possible  based  on  known  continuity  properties  across  consecutive 
segments.  These  properties  may  serve  as  constraints,  and  along  with  applicable  optimization 
methods,  may  lead  to  other  reduced  bandwidth  representations  or  to  a  useful  minimum 
bandwidth  representation. 
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List  of  Acronymns 


DFT  Discrete  Fourier  Transform 
FFT  Fast  Fourier  Transform 
FIR  Finite  Impulse  Response 
FM  Frequency  Modulation 
FSK  Frequency  Shift  Keyed 
Hz  Hertz  (cycles  per  second) 

HR  Infinite  Impulse  Response 

MP  Matrix  Pencil 
ROC  Region  of  Convergence 
S/I  Signal-to-Interference  Ratio 
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Appendix:  Sample  Pencilgram  Code 

The  following  Matlab  (The  Mathworks,  Inc.)  code  is  provided  for  informational  purposes 
only;  there  are  no  direct  or  implied  warrantees  regarding  any  aspects  of  the  code.  For  examples 
presented  in  this  report,  the  pencil  parameter  was  set  as  p  =  L/2. 


function  [X,  xr]  =  mp  gram  complex  (x,  Kc,  p,  L,  N) 


o, 

o 

o, 

o 

o, 

o 

g, 

o 

g, 

o 

g, 
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o, 

o 

o, 

o 

o, 

o 

o, 

o 

g, 

o 
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o 
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g, 

o 

g, 
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usage:  [X,  xr] 


mp  gram  complex (x, Kc, p, L, N) 


parameters :  x 
Kc 

P 

L 

N 


data  samples;  L  <=  length (x) 
estimated/specif ied  number  of  exponentials 
pencil  parameter  (choose  p  >=  Kc,  p  <=  L-Kc) 
length  of  signal  analyzed  (samples,  L  <  2000  or  so) 
number  of  discrete  frequencies  around  the  unit  circle 


returned:  X  -  unit  circle  response  (complex  valued) 
xr  -  reconstructed  estimate  of  x 

AJNoga,  December,  2008;  this  version  handles  complex-valued  inputs. 


%  unnecessary  dimensionality  change  (for  consistency  with  the  original  math) 
if (size (x, 2)  ~=  1) 
x  =  x  '  ; 

end 

%  parameter  for  ensuring  singular  pairs  are  retained  in  tact 
PCT  =  .01; 

%  check  input 
if  p  <  Kc 

error ('mp  gram  complex  -  Choose  p  >=  Kc '  )  ; 

end 

if  L-Kc  <  p 

error ('mp  gram  complex  -  Choose  L-Kc  <  p '  )  ; 

end 

%  form  the  matrix  X0;  alternatively,  can  use  the  following  line 
%  X0  =  hankel (x ( 1 : p) , x (p : L) ) ; 

X0  =  zeros (L-p, p+1 ) ; 
for  i=l:L-p; 

X0(i,l:p+1)  =  x ( i : p+i-1+1 ) ; 

end 

%  singular  value  decomposition  of  X0 
[U, Sd, V]  =  svd (X0) ; 
dS  =  diag ( Sd) ; 
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%  take  care  of  double  poles 
if (Kc+1  <=  p) 


if (dS (Kc)  ~=  0) 

if ( (dS (Kc) -dS (Kc+1) ) /dS (Kc)  <  PCT) 


Kc  =  Kc  +  1 


end 


end 

end 


o, 

o 


truncate  to  Kc  singular  values 


U0 

U1 


=  U (1 :p-l, 1 :Kc) ; 
=  U  (2  :p,  1  :Kc)  ; 


%  equivalent  to  calculating  the  eigenvalues  of  X02Kinv*Xl 

%  (row  selections  of  decomposition  rather  than  column  selections  of  data) 
lambda  =  eig (pinv (U0 ) *U1 ) ; 

%  "equation"  matrix.  A,  for  least-squares  estimation  of  amplitudes,  phases 
for  ii=l:L 

for  jj=l:Kc 

A  (ii, j j ) = lambda (jj)A(ii-l); 


end 


end 


%  least-squares  solution  matrix 

T  =  inv (A ' *A) * (A ' ) ;  %  same  as  inv (conj (A' ) *A) *conj (A' ) ; 

%  project  data  onto  "basis  functions"  to  get  the  solution,  x 
x  =  T*x; 

%  least  squares  estimates  of  phases  and  amplitudes 
r  =  abs (x) ; 
ph  =  angle (x) ; 

%  generate  X  and  signal  estimate,  xr 
X  =  zeros (N, 1 ) ; 

W  =  ( (0 :N-1) *2*pi/N) ' ; 
n  =  ( 0 : L-l )  '  ; 
xr  =  zeros (L, 1 ) ; 
z  =  exp ( j  *W) ; 
for  k  =  1 : Kc 

a  =  lambda ( k); 

alphak  =  log (abs (a)); 

Wk  =  angle (a) ; 

%  evaluate  on  the  unit  circle 
F  =  real ( z . / ( z-a)  -  .5); 
if (alphak  <  0)  %  decaying  component 

X  =  X  +  r ( k) *exp ( j  *ph ( k) )  . *F; 
else  %  growing  component 

X  =  X  -  exp (alphak* (L-l )) *exp ( j *Wk* (L-l )) *r ( k) *exp ( j *ph ( k) ). *F; 

end 

xr  =  xr  +  x ( k) *( lambda ( k) . An) ;  %  regenerate  signal 


end 
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